Plastid Phylogenomic Insights into the Inter-Tribal Relationships of Plantaginaceae

Simple Summary The current classification of Plantaginaceae divides the family into 12 tribes. However, the inter-tribal relationships of this family remain unresolved, and phylogenetic and comparative studies on Plantaginaceae plastomes are also scarce. In this study, we compared 41 plastomes representing the 11 tribes of Plantaginaceae and uncovered the high conservation of the plastomes in terms of genome structure and gene content. Phylogenetic analyses based on 68 plastid protein-coding genes (PCGs) successfully elucidated the inter-tribal relationships of Plantaginaceae, especially disputed tribes, such as Plantagineae, Digitalideae, Veroniceae, Hemiphragmeae and Sibthorpieae. PhyParts analysis showed a certain extent of conflict between gene trees and species tree, revealing the limitations of phylogenetic analysis of Plantaginaceae using single or multiple plastid DNA sequences. Collectively, our results provide some basic information on the plastomes of the Plantaginaceae taxa and some new insights into the inter-tribal relationships of Plantaginaceae, laying the groundwork for future phylogenetic and taxonomic studies. Abstract Plantaginaceae, consisting of 12 tribes, is a diverse, cosmopolitan family. To date, the inter-tribal relationships of this family have been unresolved, and the plastome structure and composition within Plantaginaceae have seldom been comprehensively investigated. In this study, we compared the plastomes from 41 Plantaginaceae species (including 6 newly sequenced samples and 35 publicly representative species) representing 11 tribes. To clarify the inter-tribal relationships of Plantaginaceae, we inferred phylogenic relationships based on the concatenated and coalescent analyses of 68 plastid protein-coding genes. PhyParts analysis was performed to assess the level of concordance and conflict among gene trees across the species tree. The results indicate that most plastomes of Plantaginaceae are largely conserved in terms of genome structure and gene content. In contrast to most previous studies, a robust phylogeny was recovered using plastome data, providing new insights for better understanding the inter-tribal relationships of Plantaginaceae. Both concatenated and coalescent phylogenies favored the sister relationship between Plantagineae and Digitalideae, as well as between Veroniceae and Hemiphragmeae. Sibthorpieae diverged into a separate branch which was sister to a clade comprising the four tribes mentioned above. Furthermore, the sister relationship between Russelieae and Cheloneae is strongly supported. The results of PhyParts showed gene tree congruence and conflict to varying degrees, but most plastid genes were uninformative for phylogenetic nodes, revealing the defects of previous studies using single or multiple plastid DNA sequences to infer the phylogeny of Plantaginaceae.

tree [28]. The objectives of our study are to (1) investigate the structural and compositional variation of plastomes within Plantaginaceae; (2) reconstruct a robust phylogeny and preliminarily elucidate the relationships of the disputed tribes Plantagineae, Veroniceae, Digitalideae, Hemiphragmeae and Sibthorpieae; (3) assess the degree of gene-tree concordance and conflict within plastid data of Plantaginaceae, revealing the limitation of phylogenetic analyses using single or multiple plastid DNA sequences.

Plant Materials, DNA Extraction, Sequencing, Assembly and Annotation
A dataset of plastomes from 43 species was used in this study, representing 11 tribes of Plantaginaceae (except Globularieae). Among them, six species (Trapella sinensis, Adenosma glutinosum, Deinostema violacea, Ellisiophyllum pinnatum, Russelia equisetiformis and Linaria buriatica) were newly sequenced in this study. All voucher specimens were deposited in the Herbarium of Kunming Institute of Botany (KUN), Chinese Academy of Sciences. Voucher information, the source of material and GenBank accession numbers are given in Table S1. Pastomes of the other 37 species (including 2 outgroup taxa) were obtained from GenBank of the National Center for Biotechnology Information (NCBI) Table S2).
Genomic DNA of each sample was extracted from the silica gel-dried leaf tissues using the modified CTAB method [29]. Paired-end libraries with an average insert size of approximately 400 bp were prepared using a TruSeq DNA Sample Prep Kit (Illumina, Inc., San Diego, CA, USA) according to the manufacturer's instructions. The libraries were sequenced on the Illumina NovaSeq platform at Personalbio (Shanghai, China). Raw data were filtered using fastP v0. 15.0 (-n 10 and -q 15) to obtain high-quality reads [30]. Clean paired-end (PE) reads were assembled using GetOrganelle v1.7.1 [31] with K-mer values of 21,45,65,85, and 105 as seen using SPAdes v3.10 [32]. The complete plastome of Bacopa monnieri (MN736955) was set as the reference for Trapella sinensis, Adenosma glutinosum, and Deinostema violacea, while Antirrhinum majus (MW877560) was used for assembly the plastomes of Linaria buriatica, Russelia equisetiformis, and Ellisiophyllum pinnatum. Bandage v0.8.1 [33] was used to visualize and filter out the assembled contigs to generate a complete circular plastome. Plastid Genome Annotator [34] was used to annotate the newly assembled plastomes and then manually check the consistency of start/stop codons against the reference genome in Geneious v10.2.3 [35]. The tRNAscanSE web service [36] was used to confirm the tRNA genes with default parameters. Furthermore, the online program OrganellarGenomeDRAW (OGDRAW) [37] was used to draw circular plastome maps. GenBank accession numbers of the six plastomes generated in this study are presented in Table S1.

Phylogenetic Analyses
Typically, organellar genomic regions are concatenated as a single locus to reconstruct phylogeny. However, a recent study uncovered notable levels of gene-tree conflict within the plastome [40], suggesting that concatenated analyses may be inappropriate for plastomes and should be performed with caution. In light of this, Gonçalves et al. [41] advocated for the use of both concatenation and multi-species coalescent methods (MSC) for inferring plastid phylogeny. Multispecies coalescent algorithms are divided into summary methods that use estimated gene trees to infer the species tree and "single-site" methods that use nucleotide alignments for species tree inference. Despite not performing full coalescent analysis, the summary method ASTRAL [27], used to estimate a species tree given a set of gene trees, has been shown to be statistically consistent under the multi-species coalescent model [42]. Therefore, we utilized both concatenated and ASTRAL analyses to infer phylogenetic relationships.
Before phylogenetic analyses, all 68 protein-coding genes (PCGs) were extracted from the plastomes of 43 taxa (Tables S1 and S2) using Geneious v10.2.3 [35]. Each proteincoding gene was separately aligned using MAFFT v7.450 [43] with the default settings and adjusted manually. Characteristics of alignments of 68 PCGs involved in the phylogenetic analyses are listed in Table S3. Scrophularia dentata and Buddleja sessilifolia (Scrophulariaceae: Lamiales) served as the outgroup according to the results of Gormley et al. [13].
For concatenated analysis, the multiple alignments of genes were concatenated into a 61,841 bp matrix (Table S3) with a preparation of a partition file generated by Phy-loSuite v1.2.2 [44]. Maximum likelihood (ML) analyses were performed by IQ-TREE v2.1.2 [45] using UFBoot2 [46] with the most suitable partition models (Table S4) found by Modelfinder [47] with 1000 replicates. Besides, the SH-aLRT test [48] was used to assess branch supports. According to the Akaike information criterion (AIC), Partition-Finder v2.1.1 [49] was used to select the best-fit partitioning schemes and models for each gene with the default values (Table S5). Bayesian inference (BI), implemented in MrBayes v3.2.7a [50], was constructed with an average deviation of split frequencies below 0.01. Approximately 1000,000 generations were conducted for the matrix, and each set was sampled every 2500 generations with a burn-in of 25%. The software FigTree 1.4.4 (http://tree.bio.ed.ac.uk/software/Figtree/ (accessed on 1 September 2022)) was used to visualize the phylogenetic tree.
For the ASTRAL analysis, 68 plastid gene trees were inferred separately using RAxML v8.2.11 [51] with the GTRGAMMA model and 100 bootstrap replicates, and low-supported branches (<10% bootstrap support) in the gene trees were collapsed to improve species-tree inference [52]. Subsequently, the collapsed gene trees were imported into ASTRAL-III [27] with the default settings, which produced an estimated topology of the species tree with branch lengths and local posterior probabilities (LPP) as branch support values.
In addition, PhyParts [28] was used to examine how individual gene trees agree/conflict with the species tree by mapping the 68 plastid gene trees onto the species tree generated by ASTRAL analysis, with bootstrap support (BS) < 70% at gene-tree branches considered uninformative. A Python script by M. Johnson (https://github.com/mossmatters/ phyloscripts/tree/master/phypartspiecharts (accessed on 20 October 2022)) was used to visualize the output of the PhyParts analyses. Pie charts on the species phylogeny show the number of gene trees that were concordant, conflicting, or uninformative concerning each node in the species tree.

Features of the Plastomes
Illumina sequencing generated 6,069,200-30,408,958 paired-end clean reads for the six samples. Among them, 302,408-5,093,971 reads were mapped to the final assembly, with the average coverage ranging from 295.979× to 4,997.812× (Table S6). All newly sequenced plastomes were assembled into a typical quadripartite structure containing a large single copy (LSC) and a small single copy (SSC) separated by two inverted repeats (IRs, including IRa and IRb) ( Figure 1).  (Table 1), which identically consist of a pair of IRs (IRa and IRb, with lengths of 21,404-38,724 bp) separated by a LSC region (77,713 bp) and a SSC region (4577-18,447 bp) ( Table 1). The total GC content ranged from 37.4% to 39.0%, with the lowest GC in Trapella sinensis and Limnophila sessiliflora and the highest GC content in Littorella uniflora (Table 1). Most plastomes contained the same 114 unique genes, including 80 protein-coding genes, 30 tRNA genes and 4 rRNA genes (counting all duplicated genes only once) ( Table 1, Table S7). Notably, the ycf15 gene is lost from the plastomes of Plantago, Littorella uniflora, Veronica polita, Veronica persica, Veronica eriogyne and Angelonia angustifolia. In addition, Littorella uniflora lacks functional copies of ndh genes (except ndhE); Angelonia angustifolia and Scoparia dulcis lack the infA gene, and a novel tRNA gene (trnL-CAG) is present in plastomes of Plantago maritima (Table S7).   Table 1). The total GC content ranged from 37.4% to 39.0%, with the lowest GC in Trapella sinensis and Limnophila sessiliflora and the highest GC content in Littorella uniflora (Table 1). Most plastomes contained the same 114 unique genes, including 80 protein-coding genes, 30 tRNA genes and 4 rRNA genes (counting all duplicated genes only once) ( Table 1, Table S7). Notably, the ycf 15 gene is lost from the plastomes of Plantago, Littorella uniflora, Veronica polita, Veronica persica, Veronica eriogyne and Angelonia angustifolia. In addition, Littorella uniflora lacks functional copies of ndh genes (except ndhE); Angelonia angustifolia and Scoparia dulcis lack the inf A gene, and a novel tRNA gene (trnL-CAG) is present in plastomes of Plantago maritima (Table S7).

Plastome Structural Variation
The IRs/LSC and IRs/SSC junctions among the plastomes of 15 diverse species representing 11 tribes in Plantaginaceae were compared to identify IR expansion/contraction. Overall, the boundaries of the LSC/IRb, IRb/SSC, SSC/IRa and IRa/LSC are relatively conservative within Plantaginaceae plastomes, with the rps19, ndhF, ycf 1 and trnH genes located in the above junctions, respectively ( Figure 2). The collinearity analysis also showed the conservation of plastome gene arrangement in most tribes ( Figure 3). In contrast, three Plantago species, Plantago fengdouensis, P. media and P. maritima, have distinct IR expansions, as well as rearrangements, in the SSC and IRs, which leads to differences in the junction and gene arrangement (Figures 2, 3 and S1). In addition, the size reduction in the LSC, SSC, and IRs of Littorella uniflora, especially the loss of ndh genes, brings about changes in the IRb/SSC junction ( Figure 2).    Figure S1 for more details).

Concordance and Conflict of the Gene Tree
PhyParts analyses revealed the numbers of consistent, conflicting, or uninformative gene trees for each node in the species tree, with a bootstrap support (BS) threshold of 70 set at nodes of each gene tree ( Figure 6). From the pie chart, although a certain number of genes exhibited concordance, the majority of the plastid genes were largely uninformative for nodes of the phylogeny. Besides, there are gene-tree conflicts for nodes of the species tree of varying degrees. Notably, some nodes with the conflicting genes numbers close to or greater than the numbers of concordant genes seemed to have relatively low support in the species tree. The node of Digitalideae and Plantagineae exhibited more gene-tree conflict, with two concordant genes and five conflicting genes, while the node of Sibthorpieae had five concordant genes and ten conflicting genes. Besides, only five genes supported the sister relationship between Cheloneae and Russelieae, while thirteen genes conflicted with the topology (Figure 6).

Concordance and Conflict of the Gene Tree
PhyParts analyses revealed the numbers of consistent, conflicting, or uninformative gene trees for each node in the species tree, with a bootstrap support (BS) threshold of 70 set at nodes of each gene tree ( Figure 6). From the pie chart, although a certain number of genes exhibited concordance, the majority of the plastid genes were largely uninformative for nodes of the phylogeny. Besides, there are gene-tree conflicts for nodes of the species tree of varying degrees. Notably, some nodes with the conflicting genes numbers close to or greater than the numbers of concordant genes seemed to have relatively low support in the species tree. The node of Digitalideae and Plantagineae exhibited more gene-tree conflict, with two concordant genes and five conflicting genes, while the node of Sibthorpieae had five concordant genes and ten conflicting genes. Besides, only five genes supported the sister relationship between Cheloneae and Russelieae, while thirteen genes conflicted with the topology (Figure 6).

Discussion
In general, the plastomes of Plantaginaceae are largely conserved in structure and gene content. Significant IR expansion and rearrangement of plastomes in Plantago were observed in three representative species including Plantago fengdouensis, P. media and P. maritima. Indeed, IR expansion and rearrangement of plastomes in Plantago have been discussed in detail by Mower et al. [17], and not discussed at length here. Furthermore, the absence of the ndh genes, ycf15 gene and infA gene found in Plantaginaceae appears to be common in angiosperms [53][54][55][56]. Gene loss can be the consequence of a sudden mutational event or the result of a slow process of accumulation of mutations during the pseudogenization that follows an initial loss-of-function mutation [57]. The loss of ndh genes only occurs in Littorella uniflora, which may be related to its amphibious lifestyle and partial reliance on Crassulacean acid metabolism (CAM) photosynthesis [17]. The ycf15 gene is lost in Plantago, Littorella uniflora, Veronica polita, V. persica, V. eriogyne and Angelonia angustifolia. Indeed, loss of the ycf15 gene has been observed in a variety of angiosperm lineages, which may have occurred independently throughout the evolution of angiosperms [54,58,59]. Angelonia angustifolia and Scoparia dulcis lack the infA gene and in many cases, the loss of the infA gene has been considered to be independently shifted to and expressed in the nucleus [55,60].

Discussion
In general, the plastomes of Plantaginaceae are largely conserved in structure and gene content. Significant IR expansion and rearrangement of plastomes in Plantago were observed in three representative species including Plantago fengdouensis, P. media and P. maritima. Indeed, IR expansion and rearrangement of plastomes in Plantago have been discussed in detail by Mower et al. [17], and not discussed at length here. Furthermore, the absence of the ndh genes, ycf 15 gene and inf A gene found in Plantaginaceae appears to be common in angiosperms [53][54][55][56]. Gene loss can be the consequence of a sudden mutational event or the result of a slow process of accumulation of mutations during the pseudogenization that follows an initial loss-of-function mutation [57]. The loss of ndh genes only occurs in Littorella uniflora, which may be related to its amphibious lifestyle and partial reliance on Crassulacean acid metabolism (CAM) photosynthesis [17]. The ycf 15 gene is lost in Plantago, Littorella uniflora, Veronica polita, V. persica, V. eriogyne and Angelonia angustifolia. Indeed, loss of the ycf 15 gene has been observed in a variety of angiosperm lineages, which may have occurred independently throughout the evolution of angiosperms [54,58,59]. Angelonia angustifolia and Scoparia dulcis lack the inf A gene and in many cases, the loss of the inf A gene has been considered to be independently shifted to and expressed in the nucleus [55,60].
Given the weak support and conflicting topology in the phylogenetic tree, previous phylogenetic analyses based on single or several different plastid DNA sequences have failed to clarify the tribe-level relationships within Plantaginaceae [7,[11][12][13], especially the phylogenetic placements of Russelieae, Cheloneae, Sibthorpieae, Digitalideae, Plantagineae, Veroniceae, and Hemiphragmeae. Indeed, the results of our concordance and conflict analysis verify these limitations. As shown in Figure 6, due to the lack of sufficient phylogenetic signal [40], most plastid genes are considered uninformative for nodes (BS < 70%) in species tree phylogeny. In addition, against a backdrop of largely uninformative genes, a few genes exhibit well-supported conflict at nodes to varying degrees, which can also correspond to the inter-tribal conflict in previous studies.
In most cases, systematic and stochastic error (presence of a non-phylogenetic signal in the data; the length of the genes) has been invoked by multiple studies to explain the potential source of conflict across the plastome [40,[61][62][63]. Additionally, biparental inheritance and heteroplasmic recombination (both intra-and interspecific) have also been proposed as the potential source of biological conflict [40,[64][65][66][67][68]. Here, we cannot conclude the causes of the well-supported conflict observed at nodes in our phylogeny based on the present limited analyses, which is also beyond the scope of this study. However, significantly, nodes with apparent conflicts require further research into the causes of conflict, particularly the possibility of biological conflict [40], which might help to understand the evolutionary history of plastomes of related taxa in Plantaginaceae. Moreover, the non-phylogenetic signals which can swamp the faint genuine phylogenetic signal presented in phylogeny also need to be appreciated [61,62].
In contrast, our phylogenomic analyses recovered the phylogenetic relationships of 11 tribes (except Globularieae) with robust support based on 68 plastid PCGs, providing new insights for better understanding the inter-tribal relationships of Plantaginaceae. On the whole, the phylogenetic positions of Gratioleae, Angelonieae, Russelieae, Cheloneae, Antirrhineae and Callitricheae are largely consistent with previous studies [7,8,[10][11][12][13][14][15] and are strongly supported in our study. Therefore, the relationship between these tribes is not too controversial. Notably, the sister relationship between Russelieae and Cheloneae was recovered with robust support (UFboot = 100; SH-aLRT = 100; PP = 1.00), which was unprecedented in previous analyses [7,8,10,12,13,15]. Despite the substantial morphological differences between Russelieae and Cheloneae, pair-flowered Cymes (PFCs) were indicated as the common characteristic present in the ancestors of both tribes [69]. Wolfe et al. [15,70] also definitely argued that the cymose inflorescence of Russelieae-Cheloneae would be a synapomorphy of the two tribes.
The placement of Sibthorpieae was elusive in the phylogenetic tree based on different sequence data [7,13]. However, our plastid phylogenomic analyses showed that Sibthorpieae is more closely related to the clade of Plantagineae, Digitalideae, Veroniceae, and Hemiphragmeae with strong support (UFboot = 100; SH-aLRT = 99.7; PP = 1.00), which was not proposed by previous studies. Sibthorpieae comprises the two small genera Sibthorpia and Ellisiophyllum [7]. Floral morphology with four stamens, a five-lobed corolla, and a five-lobed calyx in Sibthorpieae (mainly in Ellisiophyllum) might show a little affinity with some species of its sister clade, particularly Hemiphragma and possibly Wulfenia of Veroniceae [71]. In addition, the morphological features of pollen grains of Sibthorpieae appear to be transitional. The types of exine sculptures of pollen grains in Sibthorpieae are not only common in most of the species in Russelieae, Cheloneae and Antirrhineae, but are also typical for Plantagineae and Veroniceae [13,[72][73][74][75]. Chemically, although Sibthorpieae lacks iridoids and simple phenylethanoids widespread in Plantaginaceae, some of the caffeoyl phenylethanoid glycosides (GPGs) present in Sibthorpieae are similar to those found in Plantagineae, Digitalideae and Veroniceae [76]. Therefore, the placement of Sibthorpieae seems plausible in our study. Nonetheless, the new placement of Sibthorpieae should be treated with caution due to not sampling tribe Globularieae, which used to be related to a clade containing Plantagineae, Digitalideae, Veroniceae and Hemiphragmeae [7,8,11,13,77].
Plantagineae, Veroniceae, Digitaleae and Hemiphragmeae formed the Clade VI in our phylogenetic analyses. However, the inter-tribal relationships within this clade differed considerably between analyses. Both Olmstead et al. [77] and Bello et al. [8] recovered Digitalideae and Hemiphragmeae as successive sister groups to Veroniceae + Plantagineae with marginal support, whereas Albach et al. [7] found that the phylogenetic location of Digitalideae and Hemiphragmeae interchanged with moderate support. Moreover, Estes and Small [11] identified Digitalideae and Plantagineae as successive sister groups to Veroniceae + Hemiphragmeae, which was also supported by Gomley et al. [12], but the placement of Plantagineae was not well-supported in their studies. Unlike most previous studies, our analyses favor a sister-group relationship between Plantagineae and Digitalideae, and Veroniceae is closely allied to Hemiphragmeae. These affinities are also consistent with recent studies based on plastome, mitochondrial and nuclear ribosomal data [14,17].
Plantagineae was shown to be sister to Veroniceae in most molecular analyses [7,8,77,78], and their shared tetramerous flowers, seed morphology and some similar iridoid constituents between these two tribes were often discussed to show their affinities [7,78]. However, given that the pentamerous and zygomorphic flowers could be plesiomorphic in Plantaginaceae [79], the presence of tetramerous flowers in Plantagineae and Veroniceae was assumed to be an independent shift based on mixed evidence for fusion and loss of flower parts [5,71,[78][79][80]. In addition, some similar phytochemical constituents in Plantagineae and Veroniceae seem to be more widespread in Plantaginaceae [77].
Plantagineae encompasses Plantago, Littorella and Aragoa [16], and the (nearly) actinomorphic, tetramerous corolla and the four equal stamens are considered synapomorphies for this tribe [7,78]. In contrast, Digitalideae, including Digitalis and Erinus, may retain a plesiomorphic character state with a shared long-tubed zygomorphic corolla with five lobes, a deeply five-lobed calyx, and four didynamous stamens (the fifth one is evident only in early ontogeny) [7,78]. Although there are few floral morphological commonalities between Plantagineae and Digitalideae, a similar seed morphology between Erinus and Plantagineae and a five-lobed calyx shared in Digitalideae and Aragoa might show their affinity [5,78]. Furthermore, similarities in chemosystematic characters could also help to understand their close relationship. Sorbitol, a carbohydrate with a limited distribution, appears to be the sugar characteristic for Digitalideae and Plantagineae, although Erinus contains glucose [76,81]. Indeed, the sorbitol found in Digitalis is also present in all species of Plantago [82] and Aragoa [83]. In contrast, mannitol was found to be the sugar characteristic of Veroniceae [76]. In addition, although Digitalis contains unique cardenolides and lacks iridoids, Erinus preserves the pattern of iridoid biosynthesis, retaining a phytochemical arsenal similar to that of related Plantagineae [76,81].
The sister relationship between Hemiphragmeae and Veroniceae was strongly supported in our analyses, which was also recovered by Gomley et al. [12] and Estes and Small [11], but it was rarely emphasized before. Hemiphragmeae merely contains a monotypic genus Hemiphragma from the Himalayas to China, Formosa, Philippines and Celebes with a fleshy, septicidal capsule, four stamens, an actinomorphic five-lobed corolla, a five-lobed calyx, and dimorphic leaves [1,71]. By comparison, Veroniceae with worldwide distribution comprises a large genus Veronica and the wulfenioid grade including Veronicastrum, Lagotis, Wulfenia, Picrorhiza, Wulfeniopsis and Paederota [6,7,16]. Almost all genera are characterized by two stamens (except Picrorhiza, which has four stamens). Veronica sensu lato is mainly characterized by (almost) actinomorphic flowers, a short to absent corolla tube, four-lobed corolla and four-lobed calyx, whereas the majority of genera in the wulfenioid grade have long-tubed zygomorphic flowers, a two-lipped corolla (commonly four corolla lobes) and five-lobed calyx (except Lagotis, variable in corolla and calyx number, and Wulfenia with a five-lobed corolla) [6]. Obviously, the floral morphological differences between Veroniceae and Hemiphragmeae are therefore unlikely enough to relate these two tribes. Nevertheless, their sister relationship can be justified by the following lines of evidence. First, Hemiphragmeae and Veroniceae share a common pollen type with similar patterns of sculpture of exine and aperture membranes [73]. Next, despite the differences in floral symmetry, the pentamerous calyx and corolla shared in Wulfenia and Hemiphragma might show a little affinity between the two tribes. Indeed, with several characteristics such as a two-lobed stigma, pentamerous calyx and corolla, divergent anther thecae, and colporate pollen, Wulfenia is frequently considered the most primitive representative in Veroniceae [6,84,85]. Furthermore, Wulfenia, Lagotis and Veronicastrum are supposed to be most likely to form the representatives of the early diverging branch of Veroniceae in most molecular studies [4,6,7], and some of the chemical compounds present in Wulfenia and Lagotis are very similar to those found in Hemiphragma [76]. Finally, based on comprehensive analyses [4,86,87], the possible common ancestral distribution area (the Qinghai-Tibetan Plateau (QTP) or Himalayan region) of the two tribes may also provide evidence for their sister relationships. Nonetheless, further studies are necessary to validate the relationships found here.

Conclusions
The present study enriches the genomic resources of the Plantaginaceae plastome and provides basic information on the plastome structure and gene content of Plantaginaceae taxa. Although Globularieae was not included in our analyses, the reconstructed robust plastome phylogeny, including 11 tribes of Plantaginaceae, provides some new insights into the inter-tribal relationships of Plantaginaceae. In addition, to further clarify the relationships between tribes in Plantaginaceae, future studies will benefit from using more extensive taxonomic sampling of all tribes from the plastid genome and employing more data from the nuclear and mitochondrial genomes.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/biology12020263/s1, Figure S1: Corresponding rearrangement areas of plastomes from Plantago fengdouensis, P. media and P. maritima were shown, with Plantago nubicola as the reference. Linear gene maps were drawn using OGDRAW (Greiner et Table S1: list of six species newly sequenced with voucher and GenBank accession numbers; Table S2: list of Plastomes in phylogenetic analyses from GenBank; Table S3: characteristics of alignments of 68 PCGs involved in the phylogenetic analyses; Table S4: dataset partition with best-fitting models of molecular evolution (Modelfinder) for ML analyses; Table S5: dataset partition with the best-fitting model of molecular evolution (Partition Finder) for BI analyses; Table S6: summary of Illumina sequencing for six species; Table S7: list of genes identified in the plastomes of Plantaginaceae.
Author Contributions: Conceptualization, H.Y. and P.X.; data curation, P.X., L.T. and Y.L.; formal analysis, P.X., L.T. and Y.L.; writing-original draft, P.X. and L.T.; writing-review and editing, H.Y., C.L., P.X. and L.T.; supervision, H.Y. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: The sequence data generated in this study are available in GenBank of the National Center for Biotechnology Information (NCBI) under the access numbers: OQ129603-OQ129608.